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Abstract 

Among the efficient numerical metliods based on atomistic models, the quasicontinuum (QC) 
method, introduced by Tadmor, Ortiz, and Phillips (1996), has attracted growing interest in 
recent years. Originally, the QC method was developed for materials with simple crystalline 
lattice (simple crystals) and later was extended to complex lattice (Tadmor et al, 1999). In the 
present paper we formulate the QC method for complex lattices in a homogenization framework 
and perform analysis of such a method in a ID setting. We also present numerical examples 
showing that the convergence results are valid in a more general setting. 

1 Introduction 

In some applications of solid mechanics, such as modeling cracks, structural defects, or nano- 
electromechanical systems (NEMS), the classical continuum description is not suitable, and it is 
required to utilize an atomistic description of materials. However, full atomistic simulations are 
prohibitively expensive, hence there is a need for efficient numerical method. Among the effi- 
cient methods based on atomistic models, the quasicontinuum (QC) method has attracted growing 
interest in recent years |24j . 

The QC method is a multiscale method capable of coupling atomistic and continuum description 
of materials. It is intended to model an atomistic material in a continuum manner in the regions 
where deformation variations are low and use fully atomistic model only in the small neighborhood 
of defects, thus effectively reducing the degrees of freedom of the system. Originally, the QC method 
was developed for materials with simple crystalline lattice [32] and the convergence of a few variants 
of the method has been analyzed under some practical assumptions (see, e.g., [HI [25l [22l [27] ) . The 
QC method is based on the so-called Cauchy-Born rule (see, e.g., [I8l[20l[l6l[8]) which states that the 
energy of a certain volume of material can be approximated through the deformation energy density, 
which is computed for a representative atom assuming that the neighboring atoms follow a uniform 
deformation. Later, QC was extended to materials with complex lattice (a union of a number of 
simple lattice sites) [33j based on the improved Cauchy-Born rule [31] which accounts for relative 
shifts between the comprising simple lattice sites. Examples of such materials include diamond cubic 
Si, HCP metals (stacking two simple hexagonal lattices with a shift vector) like Zr, ferroelectric 
materials, salts like Sodium Chloride, and intermetallics like NiAl. Recent developments of QC 
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for complex lattices also include adaptive choice of representative cell of complex crystals |13] . It 
appears that no rigorous analysis is available so far for the complex lattice QC. 

In the present work we propose a treatment of complex crystalline materials within the frame- 
work of discrete numerical homogenization. Homogenization techniques for partial differential 
equations (PDEs) with multiscale coefficients are known to be successful for obtaining effective 
equations with coefficients properly averaged out [7]. Finite element methods based on homoge- 
nization theory have been pioneered by Babiiska [4j and have attracted growing attention these 
past few years (see [211 HH |T7] for textbooks or review papers). Following the ideas of [7], we use 
homogenization techniques to describe the coarse-graining of complex lattice. This allows to give a 
new formulation of the QC method for complex lattice (that we will sometimes call "homogenized 
QC" (HQC) method). More interestingly, we find that there is equivalence between the discretely 
homogenized QC and the complex lattice QC based on the improved Cauchy-Born rule. We can 
then use this discrete homogenization as a framework for the description of a quasicontinuum 
method for complex materials. There are several benefits in this regard. First, in this framework 
the connection to the well-developed theory of continuum homogenization and related numerical 
methods becomes more apparent. This allows us to apply the analysis techniques developed for 
continuum homogenization [21 [15] to the quasicontinuum method for complex materials. Second, 
homogenization theory can be used to upscale the atomistic model in both, time and space, which 
makes it promising for modeling and especially analyzing zero temperature and finite temperature 
motion of atomistic materials |15l \T9\ I24| . Also, homogenization can be applied to "stochastic" 
materials, atomistic counterparts of which include polymers [6j and glasses. Last, for the finite 
temperature simulations, when materials are modeled with static atoms interacting with effective 
temperature-dependent potentials, homogenization may serve as a rigorous instrument to derive 
such potentials. 

We note that the idea of applying homogenization to atomistic media has appeared in the lit- 
erature [m [m [H dHl [6]. We also note that the method considered in this paper is essentially 
equivalent to the QC for complex crystals, being put in the framework of numerical homogeniza- 
tion0 However, the rigorous discrete homogenization procedure and related numerical method 
allow us to derive error estimates for the homogenized QC method, when compared to the solution 
of discretely homogenized atomistic equations. It also allows, by a reconstruction procedure, to 
approximate the original full atomistic solution. To the best of our knowledge, such error estimates 
are new. As in many numerical homogenization techniques for PDEs, there is no need for our 
numerical approximation to derive homogenized potential before-hand, since the effective potential 
is computed on the fiy (see, e.g., |15]). Finally, we note that the error estimates are derived in 
one dimension for linear interaction, but the numerical methods itself applies to nonlinear multi- 
dimensional problems. Numerical experiments show that the derived estimates are valid in more 
general situations. 

The paper is organized as follows. After a brief presentation of a ID model problem of atomistic 
equilibrium (Section [2]) we discuss discrete homogenization (Sections [3] and |4|), and then formulate 
and analyze a macro-micro numerical method capable of capturing effective behavior of a complex 
material (Sections [5] and [6]). We also illustrate how the presented technique can be applied to 
2D crystals (Section [7|). Numerical examples illustrating the performance of our method are then 
presented (Section [8|), followed by concluding remarks (Section [9]). 

^For more details on relations of different multiscale methods for complex crystalline materials, refer to the 
companion paper ^3J. 
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1.1 Function spaces 

We consider the space of n-periodic functions on the lattice SZ {S e R, S > 0): 

U^eri^Z) = {u:6Z^R: u{Xi) = u{Xi+n) Vi G Z} . 
and the space of n-periodic sequences with zero average: 

U^{SZ) = {ueU^,,{dZ): {u)^ = 0}, 
where the discrete integration operator {•)x is defined for u G W^^^{5'L) by 



1 " 



1=1 

Likewise, we consider the tensor product space on dxL x 62^: 

u{Xi,») = u{Xi+ni,»), u{»,Yj) = u{;Yi+n2) Vi,j G Z}, 
and discrete integration operators 

"1 1 rb2 



1=1 j=i 
A bihnear form for u,v & J7pgj,((^Z) is defined by 

{U, V)x = ^^ ^(^i) ^(^i 

and for ii,^; G U^^,{6iZ) U^%iS2Z) by 

For u G Up^j.{5Z) we introduce the forward discrete derivative Du G Up^j.{SZ) 
and the r-step discrete derivative (r G Z,r / 0) D^v G Upg^{5Z) 

u{Xi+r) - u{Xi) u{Xi+r) - u{Xi) 



n 
1=1 



rii ri'z 



DrU{Xi) = 



Xi+r ~ 



In addition to differentiation operators, we define for u G Up^^{SZ), the translation operator Tu G 

^p"er(«) 

TuiXi)=u{X,+,). 
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Then the r-step translation (r G Z) can be expressed as a power of T: 

T''u{Xi) = u{Xi+r). 

The definitions of the discrete derivative and translation generalize to functions in U^^^{8i'L) 
^per(<^nZ) by considering the partial discrete derivative and translation operators, i.e., Dx,Dx,r, Tx 
applied to u{»,Yj) and DyiDy^tiTy applied to ^(Xj,*). 

The following lemma, whose proof is trivial, will be useful: 

Lemma 1 (Discrete integration by parts). Foru,v G U^^^{5'L) the following identity holds: 

{u, DrV)-^ = — (T~^DrU, ^ . 

This identity can he written in an operator form as (-D^)* = T~^Dr. 
We finally define appropriate norms for functions v G Up^^{()Z): 

/ n XV? 

II^'IIl-H = \'"\w^-''{n) = \\Dv\\Li{n), (1) 

II II II II r-i2 II II 

1.1.1 Identification in Mpg^ 

It is clear that a function u G Up^^{SZ) can be identified with a representant u = [ui]^^-^ in Mp^j,, 
where Ui = u{Xi) (the subscript per means that Ui is defined by periodic extension ii^+n = Ui for 
all indices i G Z). We can also identify functions in UJ^{SZ) with their representants in Rp^^ with 
zero mean. We will denote this vector space as M^. In this paper we will use a product space of 
C/pgj,((5Z) with different values of 5. In such a case it is important to retain the functional notation 
for u. However, when there is no confusion, we will avoid such heavy notations and simply use 
u,DrU,Tu G Mpgj. where due to identification of Uj with u{Xi), the operators are defined as 

{DrU)i = DrUi = , {Tu)i = TUi = Ui+i, (2) 

Di will simply be denoted as D. Likewise the discrete integration and bilinear form can be written 
as 



1=1 i=l 



n 

UiVi. 



The notation uv denotes the componcnt-wisc product: uv = [ujWjJj^^ , which will enable us to 
conveniently write {u,v)- = (uv),'. A scalar a G M will sometimes be identified with the vector 
a = [a]"^^. Finally, for the norms previously defined on U^{57j), we will use the following notations 
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for V e M^J^j.: 



n ^ 



!/</ 



= - 2^ L°°(n) = max Uj , = uv iq^^x, 

\ n -^-^ / l<«<n ^ ^ ^ ' 

2 II (''^) 



When it will cause no confusion, we will omit the argument n in the norms, thus writing only 
II^IIl2> Hm, etc. 

2 Problem Formulation 

The focus of the present study is on correct treatment of atomistic materials with spatially oscil- 
lating or inhomogeneous local properties. For simplicity, we will first consider the ID periodic case 
(the 2D case will be discussed in Section [Tj). 

2.1 Equations of Equilibrium 

We describe the formulation of the problem of finding an equilibrium of an atomistic material in 
the ID periodic setting. We consider the periodic boundary conditions in order to avoid difficulties 
arising from presence of the boundary of the atomistic material. Otherwise, the boundary of an 
atomistic material, unless properly treated, would contribute an additional error to the numerical 
solution, studying which is not an aim of the present work. Nevertheless, it should be noted that 
the numerical method and the algorithm proposed in the present work can be applied to Dirichlet, 
Neumann, or other boundary conditions. 

Consider a material at the microscopic scale which occupies a domain 0. We assume that 
the position of the atoms in reference configuration is given by Xi = ei € eZ Ci fi. When the 
material experiences a deformation the atom positions become Xi = Xi + Ui. We assume that the 
displacements Ui behave periodically with a period length G N, i.e., 

lij+AT = Ui — oo < i < oo. (3) 

Setting Ui = u{Xi) (see Section n.l.ip we see that u G U^^^{eL). According again to Remark ll. 1.11 
we will identify u with u G for the discussion which follows. 

We assume that the atoms Xi,Xj interact through the pairwise potential ^ij-, which depends on 
particular atoms i and j thus allowing for modeling heterogeneous materials. Due to the assumption 
of periodic displacements we have (piJ^M,j+N = ^ij- The energy of atomistic interaction of the system 
(summed for the atoms over one period) is then 

Af oo / N Af oo 



iv oo / \ iv oo / \ 

EM = E ) = E (u - ^) + ^) 

i=l j=i+l ^ ^ i=l j=i+l ^ ^ 

N oo , s N oo 

= e E E '^*'*+^ ( + — ~ ) ^ ^ E E '^*'*+'' + r^rUi) . 

i=l r=l ^ ^ i=l r=l 
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We assume that the potential (pij{z) vanishes for \z\ large enough, so that it is sufficient to 
consider at most R neighboring atoms in the interaction energy: 

N R I R \ 

i=l r=l \r=l / i 

where ■ ^per ~^ ^per ^-^'^ introduced in the following way: 

{^r{z))i = ^i,i+r{r + rzi). (4) 
The potential energy of the external force / is 

N 



1=1 

The forces on each atom are given and considered to be independent of actual atom positions 
Xi. For the problem to be well-posed, the sum of all forces per period is assumed to be zero, i.e., 
(/). = 0. 

The total potential energy of the atomistic system is then 

U{u) = Ei^t{u) + Ee^t{u). 
In these notations the problem of finding the equilibrium configuration of atoms can be written as 

For the equations ([5]) to have a unique solution, we must additionally require that the average of 
u is zero: 

(u)^ = 0. (6) 

The equilibrium equations ([5]) together with the additional condition ^ can be written in 
variational form: find u G such that 



N 
per 

H^ = 0, (7b) 



U'{u;v) = EUiu;v) + E',M = € <er (7a: 



where 



KM = -{f,v)^, 

R 

E{^^{u;v) = ^{^'^{Dru) ,Drv). , and 



r=l 

d d 

Thus, the variational form of the problem d?]) is 

R 

J2{K(.Dru),Drv)^ = {f,v)^ V^GM^e, (9a) 

r=l 

(tx), = 0. (9b) 
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As written in this form, the variational equation (j9ap resembles the nonlinear continuum equa- 
tion 

The problem ([9]) is often solved with the Newton's method. It consists in choosing the initial 
guess u^^^ and performing iterations to find it*-"-*. For that, the equations ([9]) are first linearized on 
the solution it^"): 



r=l 

+ J] (Z),«W) = yv£R^^^ (10a) 



and then solved for the next approximation u^^~^^^ until two successive iterations give close results. 
Here, according to dU and Q, is given by 

d 



Notice that we have identified here the N x N diagonal Jacobian matrix with a vector in 
and used the component-wise product between two vectors (see Section II. ip . 



2.2 Linearized Model and Nearest Neighbor Interaction 

We can linearize the problem in a neighborhood of a given displacement Ui: 

{^r{DrU))i ^ ripi^i_^_r{r + rDrUi) + r'^^'li+r{f + rDrUi)Dr{Ui - Ui). 

Hence upon defining 

ir = V'f'i,i+r{'r + rDrUi) - r^if'li^^ir + rDrUi)DrUi] , and 

i^r=[r'vh+ri^ + rDrUi)]^^^, (11) 

we can use the linearized approximation 

R 



U'{u; v)mY^ {^^ + Ip^DrU, Drv)i - if, v). 

r=l 

r=l \ r=l I 
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Here we used the formula of integration by parts (Lemma[T]) and again the component- wise product 
for ^'j. (Dru). We see that in the case of hnear interaction, the term with can be absorbed into 
the external force /, which turns the generic equilibrium equations ([7]) to 



R 



Y,{^rDrU,Drv), = {f,v), Vt; G M^^r (12a) 



r=l 



= 0. (12b) 



If we further assume that only nearest neighboring atoms interact (i.e., that R = 1), then the 
equations ()12p are further simplified to 

{tJ^Du^Dv), = {f,v), V^GM^,, (13a) 
{u)^ = 0, (13b) 

where we denote tp = i/^i (i.e., xl^ = xj^^ for r = 1). It will sometimes be convenient to use a "strong 
form" of dl]) or ([H]), i.e., find u £ M^g, such that 



D{ipDu) = Tf (14a) 
{u)^ = 0, (14b) 



which is derived using Lemma [TJ 



3 Homogenization of Atomistic Media 

We come now to the main subject of this paper, the treatment of materials with heterogeneous 
atomistic interaction as illustrated in Figures [1] and [2l 

h h kp ki kp-i kp 

Figure 1: Illustration of a ID model problem with heterogeneous interaction. 



Naive coarse graining for such models (e.g., as given by the straightforward application of the 
quasicontinuum method) fail to give the correct answer. One way to treat such problems is to 
apply the so-called Cauchy-Born rule for complex lattices [311 [33l [30l [13] . We present here another 
coarse graining strategy based on homogenization ideas. We derive below a discrete homogenization 
of the atomistic material which will be the basis for formulating and analyzing a quasicontinuum 
method for complex lattices. We note that our approach is different from the approach chosen in 
other works discussing homogenization of atomistic media [Qj [TOt [T9] , which consists in treating the 
homogenized material at the continuum level and the heterogeneities at the atomistic level (the 
idea of continuous X and discrete Y can also be seen in the proof of the main results in [16] ) . In our 
approach, the homogenized material will retain its atomistic description. In this section we derive 
the homogenized equation using asymptotic expansion. Rigorous justification of the homogenized 
limit will be given in Section 4 by means of error estimates towards the full atomistic solution. 
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Figure 2: Illustration of a 2D model problem with heterogeneous interaction. 



3.1 Asymptotic expansion 

We will assume that the local heterogeneity of the atomistic interaction is periodic with period 
pe, p G N. In order to take into account the local variation of the atomistic interaction we think of 
the displacement as depending on a fast and a slow scale u{Xi) ~ u{Xi, Xi/e). We define Xi £ eZ, 
the macro ("slow") variable, and Yi = Xi/e £ 7j, the micro ("fast") variable, and consider functions 
ti™ : eZ X Z — )• M indexed by m = 0, 1, 2 . . . As we consider periodic local interaction (with period 
pe, p G N) we assume that the functions are p-periodic in the fast variable, i.e., they satisfy 

u"'{X,,Yj+p) = u"'{X,,Yj), 

while the behavior w.r.t. Xi is similar to the previously considered 

u"'{X,+N,Yj) = u"'{X„Yj). 

Recalling the definitions of Section \T7l\ this means u"^ G C/j^j,(eZ) (g) C/pcr(Z). We then consider the 
asymptotic expansion 

u = u^{Xi,Yj) + eu\X^,Yj) + e\\Xi,Y,) + ... (15) 

In addition to the discrete derivative, translation, and integration, defined in Section 11.11 we need 
the total derivative and the total r-step derivative of a function u G U^^^{e'L) [/per(Z): 

u{Xi+i,Yi+l) - u{Xi,Yi) u{Xi+r,Yi+r) -u{Xi,Yi) 
Uu = , L>rU = . 

e er 

A simple calculation shows that the total derivative and the total r-step derivative can be expressed 
in terms of Dx,Dy,Ty, Dx^tiDy^t-, the discrete partial derivatives and translation operator defined 
in Section II. H in the following way: 

Du{XuYj) = DxTYu{Xi,Yj) + e-'DYu{Xi,Yj), (16) 
DMXi,Yj) = Dx,rT^u{Xi,Yj)+e'^DYM^i^Yj)- 
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3.2 Homogenization for Nearest Neighbor Linear Interaction 

In this subsection we will perform the asymptotic analysis for the equation of equilibrium of atom- 
istic materials. To explain our procedure, we first treat the simplest interaction model, i.e., the 
case of nearest neighbor linear interaction. Asymptotic expansion and homogenization procedure 
for more general cases will be given in the following subsection. We consider the problem (113p . 
written in functional form (ui = u{Xi)): 

{rDu,Dv)x = {f,v)x yveU^,,{eZ) (17a) 
(^>x = 0, (17b) 

with ip^ defined as follows: 

where the function il){Xi,») G [/pcr(Z), i.e., the tensor is "p-periodic" in the Y variable. We assume 
that the function ip is uniformly positive in the following sense: 

i^iXi^Yj) > > V(Xi, Yj) G eZ X Z. (18) 

We also assume that the external force / it does not depend on Y , i.e., / = f{Xi). We emphasize 
that oscillatory external forces could also be considered. The homogenized equation would then 
depend on a proper average of the external forces. For simplicity we will not consider this case. 
We now proceed as in the "classical homogenization" O [71 128] and plug the ansatz ([T5|) in ()14ap 
(we will go back and forth from the variational formulation (jl7p to the strong formulation (I14p ). 
This gives 

- (T^^Dx + e-^Dy) ( TpDxTYvP + e-^ipDyvP + e^DxTyv^ + ipDyu'' 

+ e^ipDxTyu^ + eDyv? + . . . ) = /. (19) 



Here we used the identity (fT6|) . and Lemma [T] to compute the adjoint {DxTy + e ^Dy)* = 
{Tx^DxTy^ + e-^Ty^Dy). We thus obt ain a cascade of equations and collect powers of e. 
Collect the O^e'^) terms in ([19]): 

-Dy {ipDyU^) = 

is p-periodic in Y. 

Thanks to (|18p we have Dyu^ = 0. This implies that is independent of Y and only a function 
of X, i.e., 

u\X,,Y,) = u\X,). 
We next collect the 0(e~^) terms in ([T9|) : 

- Dy {ipDyu^) = Dy{ipDxu^) (20a) 
is p-periodic in Y, (20b) 

where we have used the fact that n° does not depend on Y , which implies DxTy{ipDyvP) = and 
Tyu^ = u^. As usual in homogenization we take advantage of the separation of variables of the 
right hand side of (I20ap and we let x = xi^i'i ^j) be the solution of 

- Dy {ipDyx) = Dyip (21a) 
is p-periodic in Y. (21b) 
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In view of (jlSp . this problem has a unique solution (up to an additive constant) if and only 
if {Dyip) = (solvability condition) which indeed holds due to the periodicity assumption on 
ip. Existence and uniqueness follows from the Lax-Milgram theorem for the following variational 
problem (see Lemma[2]): find x[^i-,*) ^ U^{Z) such that 

{iPDYX,DYs)y = -{iP,DYs)y ys = s{Yj) €Ul{Z). (22) 

It is then readily seen that u^{Xi,Yj) = xi^u Y^j)Dxu^{^i) solves (f20]l . The general solution of this 
latter equation involves a constant depending on Xi determined by the condition (x(^j;»))y = 
(recall that functions in the space f7^(Z) have zero average), i.e., 

u\X„Y,) = x{Xf,Yj)Dxu\Xi) + u\X,). 

Finally, collecting the 0{e^) terms in (|19p gives 

-Dy {^PDyu^) = Dy [i^DxTyu^] + T^^Dx + Dyx)Dxu'^) + f 
is p-periodic in Y. 

The solvability condition for the existence of a solution reads 

{Dy {^DxTyu^) + T-^Dx (V(l + Dyx)Dxu^) + f)Y = 0, 
leading to the homogenized equation 

-Dxi^^'Dxu'') = Txf (23a) 
= 0' (23b) 

where we choose, as for the original problem (I17p . the periodic boundary conditions and where 

= {^{1 + Dyx))y- (24) 

Thus, we obtained the equation for the homogenized displacement with the homogenized 
discrete tensor . The homogenized tensor no longer depends on the fast variable Y and 
therefore we can apply the standard QC method to the homogenized equation (I23p . The equation 
(j23|) has to be supplemented with boundary conditions. Our choice of periodic boundary conditions 
for the displacement (see ([3])) leads to searching for G U^^-j.{e'L). 

In the simple case of nearest neighbor linear interaction, the homogenized discrete tensor il^^ = 
{ip (1 + Dyx))y '^^'^ be found analytically. Indeed, from (I21ap we see that ^(1 + Dyx) does not 
depend on Y: 

ij{Yj){l + DYX{X^;Y,)) = CiXi), (25) 

from where we find 

DyX = ^- 1- (26) 
The constant of integration C = C{Xi) can be found by averaging ()26p over Y: 

0={C/i^-l)y = C{l/ij)y-l, 
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from where we find 

C = {my\ (27) 

and tlie homogenized tensor is thus 

= (V (1 + dyx))y = + ^ " )^ = ^'^^y = ^ = • 

Thus the homogenized equations (|23p are written as 

Dx ({imy^) Dxu^ = f. 



We emphasize that this procedure and the obtained results are well-known for PDEs [7, Chap. 1]. 
3.3 Generalizations 

Below we generalize the results of the previous subsection to the cases of finite range (i.e., R > 1) 
linear (Section 13.3. ip and nonlinear (Section I3.3.2p interaction, omitting this details of technical 
nature. 

3.3.1 Finite Range Linear Interaction 

One technical difficulty in this case is that there are R different differentiation operators Dr- Then 
a straightforward generalization of the results of the previous subsection would yield (Xj , Yj ) 
depending on R discrete "macroscale" derivatives Dx,rU^ ^ r < R). This approach would also 
essentially differ from the results in continuum homogenization. Therefore, instead of the identity 

Dr = Dx,r + e-^Dy^r (29) 

the following approximate (accurate to 0(e)) identity should be used: 

Dr c^T^Dx + e'^Dy^r- (30) 

The accuracy of 0(e) is enough if one seeks to obtain only the homogenized solution vP{Xi) and 
the leading term of the correction ex{Xi; Yj)Dxu^{Xi). 

The procedure can now be performed similarly to the nearest neighbor interaction case: we 
plug the ansatz (jlSp in ()12ap . As previously, we obtain that u^{Xi,Yj) = u^{Xi) and the 0(e^^) 
terms yields 

uHX,,Yj) = x{Xf,Yj)Dxu'>{X,)+u\Xi), 

where x '^^ ^ solution of 

R R 

- ^ -Dy,r (iprDy^rX) = Dy^rlpr (31a) 
r=l r=l 

X is p-periodic in Y. (31b) 

Appropriate conditions on tpr are required to ensure that (j3ip has a unique solution. Collecting the 
0(e*^) terms using the solvability conditions for (similarly to the nearest neighbor case) yields 
the homogenized equation 

-Dx = Txf, 
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where the homogenized tensor is defined as 

R 

V'" = (1 + DY,rX))Y ■ 

r=l 

A remarkable feature of the result of homogenization of material with finite range interaction 
is that the material thus homogenized contains only the nearest neighbor interaction despite the 
original model having longer interactions. This is a consequence of our choice of D in the form (|3Up . 
If we would have chosen the exact relation (j29p . then the homogenized material would contain the 
same number of interacting atoms. 

3.3.2 Finite Range Nonlinear Interaction 

In this subsection we further generalize the results to the case of a general nonlinear material ([9]) 

R 

r=l 

= 0. 

We assume that the nonlinear tensor (<I'^)'(z)j has the form ^'^.{z; Xi, Xi/e) and set Yi = Xi/e 
as previously. In accordance with the definition ([8]) it means that the interacting potential (pij 
depends on Xi and Yj and 

^Uz;X„Yj)=r^{r + rz;Xi,Yj). 

We again proceed with the asymptotic expansion. We use the ansatz (llSp (we directly assume that 
= u^{Xi) in order to simplify the argument), the approximation (j30p in the above nonlinear 
equation and identify the power of e. This yields the homogenized equation 

-Dx [(<I>°)' {Dxu^)] = Txf, 

where 

R 

r=l 

The function xi^) = xi^^ ^i^^j) solves the parametric problem 

R 

- ^y,r [(^r)'(^ + DyM^))] = (32a) 

r=l 

X is p-periodic in Y. (32b) 

Of course, structure assumptions on {^f.)' are needed to ensure a unique solution of (p2]) . We will 
not go into details here as our analysis in Section 4 and 5 will be limited to the linear case. 
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Remark 1. It is useful to highlight one more feature of the homogenized equations, before we 
proceed with describing the numerical algorithm. In terms of atomistic interaction potential (friz) = 
ipr{z', Xi,Yj) the homogenized tensor (^^Y is 

r=l ^ 

This representation can be interpreted as the averaging of the functional derivative of the original 
energy 

R 

Y'Pr{r + r{Dxu° + Dy-rX)) 

r=l 

at the corrected solution Dxu^ + Dy^rX- This fact is important in showing the equivalence between 
QC applied to homogenized material and QC for complex lattices f33^ and will be proved in the 
companion paper 




4 Analysis of Equations 

In this section we show that the original and the homogenized problems of equilibrium of materials 
with spatially oscillating properties are well-posed and that the difference between their solutions 
is 0(e) in the appropriate norms. We limit our analysis to the case of nearest neighbor linear 
interaction in the ID periodic setting, but allow the material properties to vary. Such interaction 
corresponds to the nonlinear interaction linearized on a given non-uniform deformation. 

After defining the appropriate norms for measuring the error (Section 14. ip . we state the main 
theorems (Section 14. 2p followed by proof of technical lemmas (Section 14. 3p . 

In this section, by Co, Ci, C2, C3 we denote generic constants which may depend on c^, C^, 
C^, and p, but are independent of e. 



4.1 Preliminaries 

Let u G (eZ) be the solution of (jl7p . We assume as in the previous section that the tensor if;'' 
can be written as 

r{Xi) = ij{Xi,Xi/e) = i;{Xi,Yi), (33) 

where the function ^p(Xi,•) £ C/per(Z) (i.e., is "p-periodic" in the Y variable). This holds, for 
instance, if we linearize the interaction the original nonlinear model on the displacement (cf. 
pT]) ) that can be expressed as = u^{Xi, Xi/e). We also assume that ip'^ satisfies 

< < 'il^iXi,Yj) < \/Xi G eZ, Yj € Z, (34) 
\\DxnL^(N,p) < (35) 

Let G (eZ) be the solution of ()23p where the homogenized tensor is given by 

V'°(x,) = (i/v(x„.))^S 
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and x{Xi,») G [/per(Z) is a solution of (j2ip . It clearly follows from ()34p that is also coercive 
and bounded, i.e., 

< < < e eZ, (36) 

The existence and uniqueness of a solution of (fT7|l (fMl) . and ([2T]) follow from standard argu- 
ments. For the sake of completeness we briefly sketch the proof. 

Proposition 1. Let ip satisfy (I34p and assume {f)x = 0- Then the problems (I17p and (I24p /laue 
unique solutions u, u^ G C/^ (eZ) respectively, and the following estimates hold 

Hm<c^^\f\H-^: (37) 
<c^V|h-i- (38) 

Proof. We first notice that, thanks to the condition {f) x — 0) ^ linear form on (eZ). 

Problem ()17p can then be rewritten as follows: find ti G C/^ (eZ) such that 

{i,'Du, Dv)x = if, v)x yv G {el). 

Using (j34p we have {ipDu, Du) ^ ^ c^|'u|^i and the Lax-Milgram theorem concludes the proof. 
The proof of (138p follow the lines of the above proof using (I36p . □ 



Proposition 2. Let satisfy (p^ . T/ien ([2T]) /las a unique solution xi^i,*) G f7^(Z). Moreover, 
XGt/p^,(eZ)®C/^(Z). 

Proof. The problem (j2ip can be written as follows: find x(-'^i)») G [/^(Z) such that 

(V(X„.)Z)yx,^y5)y = -(^yV'(^i,-),s)y VsGC/^(Z). (39) 

As ip is p-periodic in the Y variable, we have {DyipiXi, •))y = and the existence and uniqueness 
of a solution (depending on Xi) can be established as in Proposition [TJ Notice that x depends on 
Xi. As the equation ([39|) remains unchanged when ^p{Xi,•) is changed to ^{Xi^N ,•), we also have 

Define now the corrector 

u%X.i) = u^iXi) + ex{X^,Xi/e)Dxu\X.i). (40) 

In the following subsection we show that {u'^ — u\^i < Cie||/||^2 (Theorem [T|) and Wu'^ — u\\^2 < 
Cae 11/11^2 (Theorem[2]). 

4.2 Main results 

We start with formulating the following two technical lemmas that will be proved in Section 14.31 
Lemma 2. Let x be the solution of (j2ip . 
(a) If dMl) holds then 

II II ^ PCi^ /.IN 

IIxI|l-(7V,p) < ^T — • (41) 
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(b) If both (IMI and §5i) hold then 

iV'^lw^i.-C^v) < — and (42) 

c 

\\Dxx\\LooiN,p)<P—- (43) 

In what follows, a function of two variables (e.g., x = xi^i^^j)) may be identified with a 
corresponding function of one variable (x = x{^ii ^i/ ^))- Whenever it may cause confusion we 
will explicitly specify the function space or the norm (i.e., ||xIIl2(;vp) or ||xIIl2(;v)). 

Lemma 3. Let vP G and x S C/j^j.(eZ) (g) U^{Z) be the solutions of (p3]) and (pTj) . respectively. 
Assume that (134^ holds and that N/p G N. Then the corrector defined in (I40p belongs to C/j^j.(eZ) 
and its average is estimated as 

\{n^)j,\ <e'^\\DxixDxu')\\L.^^^^y (44) 
Furthermore, the following estimate holds: 

a 



Theorem 1. Assume that {f)x = 0? A^/p G N, and i/iai (I34p and (135p /loW. T/ien i/iere exist 
constants Ci, C2 such that 

\u^-u\hi < eCi\\f\\L2, and (46) 

Kn^>xl < e'C^2||/||L^, (47) 



where is the corrector defined in (I40p . 

Proof. To show (06]) we need to estimate the right-hand side of ()45p : 
e-i|n^-n|^i < ^ ||Z)x (x(^^; >^m)^x^x°(Xi)) 



< ^ ||(i?xx(^^;>"*+i))^xn°(X,)||^2(^) 



+ ^ ||x(Xm;>^.+i)i?l«°(^i)|L2(;v) 
< — P — |Jl//-iH ^ — (-'0\\j\\l2 

where we used (|33|) . ([38]) . and (03]) to estimate DxX^ Dxu^ , and X; respectively, and also ([8 
(Lemma [H]) to estimate through ||/||2,2. Notice that we used the estimate 

\\Dlu^{X,)\\^,^^^<Co\\f\\L2, 
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which can be obtained with the help of ()36p : 

< ||^°I)V||^, = \\D {i^^Du'') - p^O)pnO)||^, , 

by estimating the terms \\D [tp^ DvP)\\^^, ||L>'0°I|l°°, and Hl^u^ll^a using (pHa]) . (|i2D . and ([3] 
respectively. 

To show ()^7|) we need to estimate the right-hand side of 

P '^'ip -1 II i^ii , P pC^p 



□ 



Theorem 2. Assuming the hypotheses of Theorem\^ there exists a constant C3 such that 

k°-n||^,<C36 11/11^., 



Proof. Using Theorem [T] yields: 



1 



<.fl.i'll/ll„-.+.^ll/lb 



□ 



4.3 Proof of Technical Lemmas 



Proof of lemma\^ In a straightforward, but very tedious calculation, one can derive, using ()26p . 
(f27ll . and (|28l) . the exact representation 



\^ p+l-2(j-/3) 1 / g{Yj 



where 5 G C/^ is defined as g{Yj) = 2ii — j for 1 < j < p. Hence (HT]l holds: 
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To show (1421) notice that 



1 ^ 

Dx^lj\X.i) = Dx 



«i "7 / Y 

= /(X.)/(AV.)( „°"t!^''' ^ ) ' 

and hence if we additionahy assume (|35p then 

1 



\Dxij'^{Xi)\ < ip\Xi)ip\Xi+i) 



HXi+i,») 



To show (H3I1 notice that 

1 



I)xx(^*, Y,) = Dxi^\X^) { g{Yj 



\ '4>{Xi,')/Y 

+ ^lj\Xi+^)lg{Y,-.)Dx- ^ 



,0/v \/ IV \ Dxtp{Xi,> 



HXi,»)'4^{Xi+i,») I Y ' 

and use ([3l|) . ([35]) . and |5f| < | to estimate 

\Dxx{X.,Y,)\< ^l,\Xi)^l^\X,^,)^ 1—^ -\ f(-7^ 

P^'-i! / 1 



2 / y 



— + 1^— =P • 



□ 



Proof of lemma\^ Under the condition N/p G N it immediately fohows from (I40p that u'^{XiJ^i^) 
hence u'' G U^^^iel). 

Denote v = xiXi,Yj)Dxu^{Xi), so that u"" = u^{Xi) + e7;(Xi, Xj/e). Since (x)y = 0, 

b/2j TV Lp/2J 

i=l j=-[p/2l+l ^ i=l j=-[p/2l+l 
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Hence express 

{v{Xi,Xi/e))^ = {v{Xi,X,/e) - v) 



. N Lp/2J 

JfY. E {v{X,,X,/e)-v{X, + ej,X,/e)) 



and estimate the terms in the parenthesis (for — \p/2\ + 1 < j < [p/2j): 
v{Xi,Xi/e)-v{Xi + ej,Xi/e) 



-J2Dxv{Xi + ek,Xi/e) j >0 

k=0 

J2 Dxv{Xi + ek,X,/€) j<0 

k=j 

j = 0, 



Lp/2J 

\v{Xi,Xi/e)-viX, + ej,Xi/e)\<e \Dxv{X, + ke,X,/e)\. 

k=-[p/2]+l 

Thus, 

N Ip/2] Ip/2] 

\{v{X„X,/e))j,\< e — Y^ Y \DxviX, + ke,X,/e)\ 

i=i j=-rp/2i+ifc=-rp/2i+i 

N Ip/2\ 

E \DxviX, + ke,X,/e)\ 

^ i=l k=-\p/2']+l 
. N b/2j 

P i=l k=-\p/2]+l 

substituting which into the definition of u'^ finishes the proof of (j44p : 

2 

\{u'^)^\ < \{u\X,))^\+e\{v{X„X,/e))^\<0+'-^{\Dxv\)^y. 



To show (|45l) we first compute, using ([25]), (127]), and (|28]) . 

i^'Dxu" = i^'iDxTy + e^^Dy) (u° + e^I^xn") 

= + i^'DvxDxu'^ + ei^'DxTyixDxu^) 

= + Dyx)Dxu'^ + erDxTrixDxu'^) 
= i^'^Dxu^ + erDxTY{xDxu^), 
Hence compute, using ()14ap . ()23ap . and the fact that Dx {u^)x ~ 

-Dxi^Dxiu'' - {u'')x - u)) = -eDx DxTy{xDxu°)) • 
Treating this as an equation for (n"^ — {u^)x — u) ^ (e^)) upon using proposition [1] one gets 

\u'-u\h^ < —\eDx {iP'DxTYixDxn"^)) b-i < ^-^\\DxTy{xDxu'')\\l2. 

Cijj Cijj 

□ 
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5 Homogenized QC for Complex Lattices 

We formulate the homogenized quasicontinuum method (HQC) — the QC method for complex 
crystalline materials — in a framework of homogenization. We introduce HQC in the general 
ID periodic case, i.e., with finite range nonlinear interaction (see ([9])), and possibly oscillating 
external force f = f"- For the case of materials with known periodic structure (i.e., crystalline 
materials), the HQC method will be equivalent to applying QC to the homogenized equations. 
We mention however two advantages of HQC. First, the method is based on the original equation 
describing the spatially oscillating material and does not rely on effective (homogenized) equation 
derived beforehand. Such strategies have proved successful in the continuum elasticity where many 
macro to micro methods averaging the effective equations on the fly have been derived (see the 
review paper [21j and the references therein). Second, we think that the HQC can be applied to 
non-crystalline materials and to time-dependent zero- or even finite-temperature problems. 

5.1 HQC Method 

Consider the problem of finding an equilibrium of an atomistic material in the general ID periodic 
case (i.e., with finite range nonlinear interaction) ([9]). The method will be presented using macro- 
to- micro framework as used in some numerical homogenization procedures [21 [T5| [23| [2T| IM] . 

5.1.1 Macroscopic afflne deformation 

Let 

X ■= {Xi =ie, i = l,...N}, TV G N, 

be the reference lattice for the problem ([9]). In the set of indices 1 < i < we choose K values 
{K < N) ii < . . . < ix and compose a macroscopic lattice 

XH.= {Xi^; i = l,...K}, 

defining the macroscopic partition of the interval 

T={5, = [X,,,Xi,^J; k=l,...K}. 

Here we fix zi = 1 for convenience (we can do so without loss of generality due to translation 
invariance) and define ik+i = iV + 1 in accordance with the periodic extension, define H}^ = 
e{ik+i — ik) (for k = 1,. . . ,K) the length of the S^, and H = max^Z/fc. We define the space of 
piecewise affine (discrete) deformations by 

U^er = {n''€ [/pl(6Z) : 

u^imsk = f^^^'f ^^(^.) + k = i...,K}, 

- ^ik ^ik+i - ^ 

and 
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5.1.2 Sampling Domains 

Inside each macroscopic interval we dioose a representative position X^^^ 

pling domain 



Sr = {Xr.X-^<X,<Xll^+pe}, Il^^ = {^eN■.Xll^<X,<Xll^+pe}, 
and define tlie operator of averaging over the sampling domain 

The sampling domain should be chosen closer to the center of the interval — ^ ''^^ if the ma- 
terial's properties vary within the interval (more precisely, if the interaction potentials for different 
groups of p adjacent atoms are different), as we will see in Theorem [6j More sampling domains per 
macro interval may be considered for higher-order macro element space . 

5.1.3 Energy and Macro Nonlinear Form 

Define the energy of the HQC method 

SkeT r=l 

where TZk {u^), defined by ()49|) . is the microfunction constrained by in the sampling domain 
S^i^"^, and ^l{z){Xi) = ipl{r + rz{Xi)) is the energy of interaction of atoms i and i + r (i.e., ipi^i+ri^) 
in the notations of Section [21 cf. (jl])). 

The functional derivative of the above energy reads 

R 

5feer r=l 

where TZ'j^{u^ ; ) is the derivative of the reconstruction TZk{u^), and (<I>,^.)'(z) = r-^(pj.(r + rz) is 
defined in accordance with 

5.1.4 Microproblem 

Given a function E Up^^, IZk {u^) is a function defined on S^^j^^ such that TZk {u^) — u^ G U^{elj) 
and 

R 

((^r)' [DrTlk (n^)) , ^rS>^^g5.op =0 Vs G U^ieZ). (49) 

r=l * 

Remark 2. When modeling essentially nonlinear phenomena (e.g., martensite-austenite phase 
transformation), one should require that the micro structure corresponds to a stable equilibrium. 
That is, one should require, in addition to (149 p . that w = IZk {'^^) ~ ^ U^{e'L) is a local 
minimum of J2r=i {^r i^riu^ + w))) ^^^^^cp JM P- 238] 
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Remark 3. In the case of linear interaction, the reconstruction TZ^ is a linear function and hence 
TZ'j^ {u^;v^^ = lZk{v^), which makes the derivative of the HQC energy ()48p take the form 

R 

5fcer r=l 

Remark 4. The functional derivative of the HQC energy (j48p can equivalently be written as 

R 

{E^^^y{u";v^) = ^ H,Y,{i'^rnDMu^)),DrV^)^^^^.^.,, (51) 

Sk^r r=l 

by noting that DrTl'^ {u^;v^) = DrV^ + (A' 7^^ {u";v^) - Drv") and that 
R 

{my [DrTZu (n^)) , [DrK - Drv"))^^^,.., = 0, 

r=l * 

in view of ()49p . Here we used the fact that Tl'p.{u^ ; ) — G C/^(eZ) which follows from taking 
the derivative oflZk{u^) — G U^{eZ). 

5.1.5 Reconstruction 

The function TZk [u^] is defined on Sj^^. We define an extension of TZk [u^) on Sk by periodic 
extension outside -S*^^^: 

n^'^(X,) = n^+7^,,(^x^)(Xi), (52) 

j-^ mod p \ , and (• mod p) is an integer value modulo p. By extending 

the function IZk {u^) in each Sk G T, one obtains a function defined on X which we denote by 

5.1.6 Variational Problem 

We define the homogenized quasicontinuum approximation as the solution G C/^ of 
where 

If the external force is not oscillating, it could instead be computed for a single representative atom. 
Well-posedness of (j53p will be discussed in Section [6] for the nearest neighbor linear interaction. 

5.2 HQC Algorithm 

The problem ()53p is nonlinear and its practical implementation is usually done by the Newton's 
method. We briefly sketch below an algorithm for solving (|53p . 
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5.2.1 Second Derivative of the Energy 

For the Newton's method we need to compute the second derivative of the energy (from (|5ip ): 

R 



Sfcer r=l 

which, applying the similar arguments as in Remark [H can be written in a symmetric form 

(55) 



Sk€T r=l 



5.2.2 Newton's Iterations for the Macroproblem 

The algorithm based on the Newton's method consist in choosing the initial guess u^^^'^ G and 
performing iterations 

until u^'^'^^^'^ becomes close to li^")'^ in a chosen norm. 

To solve the linear system (I56p for u^'^'^^^'^ — u^"^^'^ £ [/^, we choose a nodal basis wjf (1 < A; < 
K) of Up^j.. One way to satisfy the condition {u^) ^ — would be to perform all the computations 
with one basis function eliminated (e.g., to consider for 2 < k < K), and post-process the final 
solution as — (n^)^. 

The stiffness matrix of the system (I56p will thus be 



and the load vector will be 

As given by the formula ()55p we need to compute the solution of microproblem TZk (u(")'^) on each 
sampling domain S^f^^ as well as its derivative TZ'j^{u^"'^'^ ;wf^). 

5.2.3 Solution of the Microproblem 



The microproblem (I49p can also be solved with Newton's method. For that, one needs to choose 
an initial guess u^"''^^ to lZk{u^^^'^), for instance n*^"''^)(Xj) := u^'^^'^{Xi) and solve 

B. 

(($^)' (Z),n("''^)) + ($^)" (Z),n("''^)) Dr (n^^'^^+i) - n^"''^))) (Xi) = 

yx, G si^'p, 

with respect to u^"-''^^^^ constrained by u^"-''^^^^ — g C/^(eZ), until the difference between 

^{n,i/+i) g^j^j y^(n,v) -g gjj^g^^^ chosen norm. 
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After that, we can compute = 7^'^(n(")'^; tw/^) by solving 
R 

Dr [m" (l?rn("''^)) DrW,,i) (X,) = 0, VX, G (57) 

r=l 

constrained by w^^i — G f7^(eZ). Notice that the derivative of the basis functions Dj.w^ on the 
interval Sk can either be (in which case Wk^i equals zero identically), or it-^. It implies that we 
essentially need to solve the problem (I57p limited number of times (once in the ID case, or between 
d and d + \ in M'^, depending on implementation). 

Also observe that when computing TZ'^ (u(")'^; wf^), we need to invert the same linear operator 

R 

^ -Dr ((^r)" {DrU^"''^^) Dr •) as in the final Newton's iteration, which allows for some additional 

r=l 

optimization. 

5.3 Possible Modifications of the Algorithm 

First, notice that when solving for could linearize the problem on the previous iteration 

^(").^^. In that case we would have only the linear cell problems and thus we would need only outer 
Newton's iteration, but it would be required to keep the values of the micro-solution TZk{u^"'^'^) 
from the previous iteration. Moreover, even in a practical implementation of the above algorithm 
it may be required to keep the values of the micro-solution: one needs these values to initialize the 
inner Newton iterations |33] . 

Another modification could be to compute the contribution of the external force in (I54jl for 
a single atom in the case of no oscillations in f^. 

In the case of linear interaction, the algorithm becomes simpler: one does not need to do Newton 
iterations. Nevertheless, even if the algorithm in subsection 15.21 is applied to the linear problem, 
the Newton's method would converge in just one iteration. 

6 Convergence of HQC 

In this section we study convergence of the HQC method introduced in ([53]) . We analyze the method 
for linear problems and nearest neighbor interaction ()17p . We treat the external force f{Xi) in an 
exact manner. We furthermore make a slight modification to the HQC method: we assume (|33p 
and collocate the tensor tp in ()50p and (149p in the slow variable at ^jcoii in each sampling domain 
5^*^^. That is, we solve 

SfeGT 

where 

i;l,^^{Xi) :=i;{Xicou,Xi/e) VX^ G 5^"^, and (59) 

k 

{rcouDxTZk (n^) , Dxs)^^^^ro, = Vs G C/|(eZ) . (60) 
For an illustration of a collocation point X^coii and a sampling domain S^f^^ refer to Figure [3l 
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.rep , 



Figure 3: Illustration of an element Sk = [Xi^, Xi^^^-^) {ik and i^+i are indices of nodal atoms), 
sampling domain with atom indices i^^ through i^"^ +p—l (circled atoms) and a collocation atom 
i™'' (double-circled) . 



In order to obtain the second order convergence in the L^-norm, we will assume two additional 
conditions: 

X-i^ -\- X'ij^.-^—i 

X-coii < eCcoU G T, and (61) 

<c;jo. (62) 

The condition (|6ip states that the collocation point is at most 0(e) away from the center of each 

Sk- 
in this section, the constants C4, C5, etc., denote generic constants which may depend on c^, 
C^, C^, CcoU, C^o, and but are independent of e and H . 

6.1 Main Results 

Before proceeding with analysis (Section 16. 2p . we summarize the main convergence results. The 
results are formulated in terms of Cmod) the so-called modeling error, which is defined in (j7ip . 

Theorem 3. Let u^,u^ be the solutions of problems (\23\\ and ()58|) . respectively, and assume that 
(j34p and (j35p hold. Then there exist constants C4 and C5 such that 

< C4/7||/||l2, and (63) 

Ih^ - u^^, < C5H^f\\L2 + ||e„,od||L2. (64) 

Theorem 4. Under the assumptions of Theorem let u be the solution of the original problem 
p7p. Then there exist constants Cq and Cj such that 



- u\ 



L2 



Theorem 5. Under the assumptions of Theorem\^ let u be the solution of the original problem (jl7p 
and u^''^ be the reconstruction (j52p of the solution of (jSSp . where the reconstruction operator 
T^k{u^) is defined by (I60p . Then there exist constants Cg and Cg such that 



< CsH\\f\\L2, and (65) 



.H,c „,|| / tt2 



U ' — U 



^2 <C9F^||/|U2 + ||e,,od||L2. (66) 



The modeling error e^iod, defined in (|7ip . reflects the fact that we introduce some error when 
neglecting the values of the tensor ip^{Xi) everywhere outside the sampling domains S^j^^. The 
modeling error is estimated in the following theorem. 
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Theorem 6. Under the assumptions of Theorem\^ 
(a) There exists a constant C\o such that 

||emod||L2 < lemodlni < Cio-f^||/||/f-i ■ (67) 

(h ) If additionally (I6ip and ([62]) hold then there exist constants Cu and C12 such that 

||emod||L2 < -^^^ lemodli^i < {CuH'^ + Ci2e)||/||i/-i- 

(c) If the tensor 'ip{Xi,Yi) in ()59p does not depend on Xi then Cmod = 0. 
6.2 Error analysis 

We start our analysis with the following lemma which asserts that the results stated in Section 16.11 
can be reformulated in terms of the standard QC method applied to the homogenized equations 
(j23p . Recall the definition of the homogenized tensor 



lP''{Xi) = {lP{•){l + DYX{X^,•)))Y 

and define the collocated homogenized tensor V'coul^*) — ^''^(^2=011) for Xi G S^/^"^. 

Lemma 4. Under the assumptions of Theorem\^ the reconstruction u^''^ (|52p can be written as 



k 



+ ex{Xicoii,Xi/e)Du^. 

k 



Proof. Fix the element Sk and notice that the reconstruction defined by (j60p can be written as 
T^k {u^) = + w, where w satisfies 

Notice that ^^^^ and Dxu^ are constant inside each sampling domain 5^'^^. Upon substitution 
w{Xi) = w{Xi/e) = wiYi) this equation takes the form 

(V'^oii^y*, Dys)y = -Dxu" (V'^^u, DYs)y Vs G U^iZ), 

and its solution can be written as w{Xi) = w{Xi/e) = eDxu^ xi^i'^^ih Xi/(-)t cf. (p2]) . Finally, 

k 

noticing that periodically extending w in 

TZk (u^) =u" + w = u" + eDxu^xiXicou,X^/e) (69) 

k 

yields exactly n^'"^ concludes the proof of (f68|) . □ 

Lemma 5. Under the assumptions of Theorem\^ the problem (I58p is equivalent to the following 
problem: find G such that 

{^PUDu'', Dv"")^ = (/, v"")^ Vt;^ G (70) 
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Proof. We continue the argument of the previous lemma: we first fix the element Sk and differentiate 
([69D to get 

DxTZk (n^) = Dxu^ + eDxu^ DyxiXicou, Xi/e) = Dxu" (l + £»yx(^icoii, y,)) • 
Hence we compute 

(V^„ll^x7^fc(^x^)>^^g<,.cp = z)x^^^(V'coii(i + ^yx(^.|on,y.)))^ 

Finally, the following computation then shows that the left-hand side of (|58|) is equal to that of 

m- 

(^E^QCy^^^H.^H^ = J2 Hk {rconDxnk{u''),Dxv'')^^^^.., 
SkeT 

= Hk{ijl,nDxnk{u''))^^^sr.,Dxv'' 
SkeT 

SkeT 

where we used Remark H] in the first step of this derivation, and omitted the argument (Xj) or 
(X-cou) of the functions Dxu^ , Dxv^ , and ^J? since they are constant on each interval Sk G T. 

k 

Thus, (j58|) and (jTOj) are equivalent. □ 



Define the modeling error as 

emod = u^-u^, (71) 
where is the solution of the following problem: 

{^j'Du^,Dv^)^ = {f,v^)^ ^v^eU^. (72) 



Proposition 3. Solutions , of the discretized problems (|70p and (|72p exist, are unique, and 
satisfy the following estimates: 

W^'lm < c^VU-i, l^^ki < c^VIh-^- (73) 

Proof. The statement can be proved similarly to Proposition [H by noticing that c-^ < tp^ < and 
Cip < ^coU — ^V" ^^"^ substituting the original space U^{eZ) with the discretized space U^. □ 

We can now prove Theorem [6l 

Proof of TheoremlSi Part (c) of the theorem is trivial: if -0 = tl^iYi) in (I59p . then ip^{Xi) is constant, 
hence V'coU coincides with Tp^, and therefore = . 

The bound in part (a) is based on the following estimate: 



c^(Z)(n^-u^),Dt;^>^ < |(^,%i?(u^-n^),I)^;^> 



< IIV'°-V'c%llL-c;^1/lk-i b^'ki, (74) 
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where the last estimate follows from (j73p . The difference — V'coU ™ can be estimated as 
follows: for Xi £ Sk 

hence HV'" - CiilU- < ^77<^^- Taking now supremum over = 1 concludes the proof of 

part (a). 

To show (b), first observe that in (|7U|) . Du^ and Dv^ are constant on any interval G T. 
Therefore V'coU '^^^ changed to any other tensor with the same average over Sk- Hence define 

:= + - (^°>x.e5,- Since (Cu)^^,^^ = ^x^^s,' ^^^^^^^^ ® 

coincides with the solution of 

Hence we can use the arguments of (j7i|) to estimate: 

c^lu"" -u'^Ihi < ||V'°-V'°l|L-c/||/lb-i, (75) 

where 

/ max(j™",i) — 1 

= /(Xi-X.coii)D^°(X.coii) + e Yl \Xi - Xj+i\D'^^\Xj 

\ j=min(i|°",i) 

=: {Qi + Q2)x,eSk ■ 
It is straightforward to show that Qi averages up to 



XiGSk 



and can be effectively estimated using (16ip . The second terms can be estimated as 

max(i^°",i) — 1 

IQ2I < Yl \Xi - Xj+,\ < \\D^^\Xj)\\l^^HI 

j=min(i™",j) 

Thus, combining these estimates, one gets 

- V'^^m)! < eC7eoll|I)/(X.con)| + I|l- • 



Taking maximum over all Xm, 

^"'^ ■ 2- 

and substituting — into ([75]) yields the desired result. □ 
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In view of Lemma O we will turn to analysis of the problem ()72p . We next introduce the 
(homogenized) energy norm in the space (e^): 

|2 



W 



,0 = {i)^Dw,Dw)^. 



Obviously, under the assumption due to the estimate ([36]) . the energy norm is equivalent to 
the ff^-norm: 



1 < 



,0 <C, 



^1 



|2 

1^1- 



(76) 



Lemma 6. Under the assumptions of TheoremUl let G U^{e'L) he the solution to the exact 
homogenized equations (j23|) . and e be the solution to the QC equations ([72]) . Then 
the best approximation to the exact solution u^ in the energy norm, i.e., 



IS 



^1^0 < Ih^ 



1^0 



(77) 



"iv^ e C/f , 



□ 



Proof. The result follows from 

{^^D{u'' -u''),Dv'')^ 

which states that is the orthogonal projection (in the energy norm) of onto U^. 

Following the standard procedure for the analysis of finite element methods (FEM) we will 
estimate — ///n'^||^0) where Ihu^ is the nodal interpolant of u^. This interpolant is defined 
for every function v G C/j^j.(eZ) in the following way. For a partition of X define a function 
Ihv G U^^^{eL) such that 



Then set 



Ihv{X,^) = v{X,^), k = l,...K. 
Ihv = Ihv - {inv , 



X 



Thus, for V G U^^-j.{e'L) we have Ihv G U^. 
Lemma 7. The following estimate holds: 

\v - Ihv\jji < 



2^/3 



where H = max Hi, and Hh = eiXi,,^ — Xi,). 
l<k<K + 



(78) 



Proof. By noting that \v — Invlfji 



\v - Ihv\hi 



V — Imv] , the result follows from 

K «fc+i-l 

J2 \D {v{X,) - iHviX,)) \' 

k=l i=ik 



< ^Y.-i Yl \d' (viXi) - iHv{Xi 

k=l i=ik 
2 K «fc+i-l 



6 



6 



\v\ 



) 
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where we used the discrete Poincare inequahty ([85]) in the first step (notice that 
D lv{Xi) - lHv{Xi) \ = from ([78])), and the fact that D'^Ihv = in the second step. □ 



Proof of TheoremlM The estimate (|63p (convergence in the H^-novm) follows from (j76p . (177p . 
Lemma [3 and (j67p . To prove (j64p (convergence in the L^-norm) we use the standard duality 
arguments. Consider 

Then 



= (nO - ^i^) , Z)^i;0>^ = (^/° - ^i^) , (ti;° - w"") > 

hence 

- U^Wl^ < - tt^||L2 + \\U^ - U^\\l2 < C^CiH^\\f\\L2 + ||e^od||L2- 

□ 

Proof of Theorem Follows immediately from Theorems [2] and [3l □ 

Lemma 8. Under the assumptions of Theorem\^ let u^''^ he the reconstruction ()52p and u'^ he the 
corrector defined in PU]) . Then there exist constants C13 and C14 such that 



H,c „,H,c 



n^'%i< Ci3H\\f\\H-i, and (79) 



IL2 

where 



^H,c_^H,c^^^^ < Ci4 6^11/11^-1, (80) 
H,c „.H , „ ./ /,\ 7-i„,H 



and X is defined as a solution to (j2ip or ()39p J 
Proof. Notice that 



e(x - XcoiO-Dn^, 



where collocated x is defined as Xco\\{X j,Yi) = x(^jcoii, Yi) for G S'^, and can be estimated using 
iSD as 

c 

\x{Xj,Y^i) - Xco\\{X,,Y.i)\ < iX.coii - Xj\ \\Dx\\l^ < Hp^. (81) 



Then, using (|8T|) and (|73]) we obtain (|80]) : 

from where (j79p follows directly by applying the inverse discrete Poincare inequality (j86p . □ 
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Lemma 9. Under the assumptions of Theorem\^ let be the corrector defined in ()40p . Then 
there exist constants C15 and Cig such that 



uH''^-u^\^^<Ci,H\\f\\L2, and 
r,H,c „,c II ^ r~< i-r2 1 



\u'''''-u'\\^,<CieH'\\fU2 
Proof. We express 

n^.- -u'= [u^ - + [exD (n^ -u')], 
and estimate the second term of the right-hand as 



\u — u 



using the inverse discrete Poincare inequahty ()86p and the estimate (j4ip . The result follows then 
from Theorem [3l □ 

Proof of Theorem O The inequalities (j65p and (|66p follow from Lemmas [H [U and [T^ the fact that 
eH < H^, and Theorem O □ 



7 Example of Application to a 2D Lattice 

In this section, we consider a simple 2D model to illustrate how the proposed approach can be 
applied to 2D materials. 



7.1 Notations 

All the coordinates and atom indices will be vectors with two components, for instance Xi = 
{Xi^i, Xi^2)i i = (^15*2)- The unit vectors in our 2D space will be denoted as ei = (1,0) and 
62 = (0, 1). The length of a 2D vector v is denoted as \v\ = [vf + V2Y^'^, the scalar product of two 
vectors v and w is denoted as v ■ w = viwi + V2W2- 

We define the inequalities for 2D vectors in the following way: u < w if, by definition, ui < vi 
and U2 < V2 (likewise for relations >, <, and >). Thus, (1,1) <i<N means 1 < ii < A^i and 

N 

1 < ^2 < A'^2- We will also use the associated notations for the sums, for instance ^ • . 

i=(i,i) 



7.2 Equations of Equilibrium 

Consider a square lattice with the reference configuration of atoms given by 

Xi = ei ((1,1) < i < iV). 

The position of the atoms Xi and displacements Ui are related through Xi = Xi + m. We consider 
the system with A^-periodic (A^ = (A^i, A2)) conditions Ui^N = Ui. The space of such A-periodic 
vector- valued sequences is denoted as C/^r- 
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Consider the following operators on U^^^: 



DaUi = [a = 1, 2), DrUi - 



H — TT 
e r 



the averaging in i: 

1 ^ 

and the scalar product 

1 ^ 

{u,v)^ = {u.v).^ = —- ^i-^i- 

Consider the linear interaction of atoms defined by a set of neighbors IZ so that the functional 
derivative of the interaction energy is 



where V'r.i defines interaction between atoms i and i + r. Such linear interaction corresponds to a 
spring model with zero equilibrium length (cf. [20] for the discussion on the nonlinear model with 
ideal springs but with nonzero spring equilibrium length). The derivative of the external potential 
energy is 

Thus, the equilibrium equation has the form 

Y,{-4^rDrU,Drv), = {f,v)^ G C/p^,. (82) 

7.3 Homogenization 

Homogenization of equations (j82p follows Section [3l By analogy with the ID case, we will use the 
term "vector-valued function" (or, in short, "function") for u = u{Xi,Yj) rather than the term 
"vector field". 

7.3.1 Fast and Slow Variables 

We first define the fast variable Yi = Xi/e, the differentiation operators 

JJx,rV{Xi,Yj) - — — , Dx^ - L>x,e^ 

\^i+r — 



where r S Z^, r 7^ 0, a = 1, 2, the translation operators in Y: 

TY,rV{Xi,Y,) = v{X„Y,+r), TY^v{Xi,Y,) = viX„Y,+,J, 
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and averaging and scalar products: 



N 

= ITO u{X„Yj), {u,v)x = {u-v)x, 
i=(i,i) 

^^^Y = ^ E u{X„Yj), {u,v)y = {u-v)y , 
i=(i,i) 

{u)xY = {{u)y)x ' ^)XY = ■ v)xY ■ 



The space of p-periodic functions {p = {pi,P2)) is denoted as Uper and hence the functions of X 



and Y belong to C/i^j. (E> Upe: 



7.3.2 Solution Representation 

Assume, as before, Ui = u{Xi,Yi) and 'ipr^i = ipriXijYi). Then Dr = Dx,rTY,r + £~^DY,r for such 
functions u and ipr- Same as in ID case, use the following 0(e) identity: 



2 

Dr ~ ^Dx, + ^Dx, + e-'DY,r = V "^Dx^ + e-'DY, 

\v\ \v\ 1 7^ 



a=l 



Then the equation ()82p takes the form 



2 2 



E ( 2: R^^^ + ^''Dy^r) n, ( 5: + e-'DY,r) v) ={f, V)^y 

reTZ \ a=l ' ' a=l ' ' / XY 



We substitute u{Xi,Yj) = u°{Xi,Yj) + eu^{Xi,Yj) + e2'u2(Xi, Yj) + O(e^) into ([83]) and collect the 
respective powers of e. As before, by collecting the 0(e~^) and the 0(e~^) terms we obtain that 

2 

Y,) = n°(Xj) and ^^(Xi,^,) = Yl XaiXi;Yj)Dx^u^{Xi) + u^{Xi), where the matrix- valued 

a=l 

functions xi and X2 are defined as solutions to 

^ {iJrDY,rXaep, DY,rS)Y = " E ( ^^'77^^' ^^'^^ ) ^Per' = 1' 2- 

r.=l r=l \ \^\ I Y 

Collecting the 0{^) terms yields 

2 2 2 
«=1,9=1 /3=1 



where the homogenized tensors ijj^g are defined as 



a/3 



y 



with / denoting a 2 x 2 identity matrix. The homogenized tensors -i/^J^^ are related to the fourth-order 
stiffness tensor in linear elasticity theory [ll [26l [29] . 
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7.3.3 Example of Application of Homogenization 

To illustrate how the 2D discrete homogenization works, we apply it to the following model problem. 
The set of neighbors is defined by 7^ = {(1,0), (1, 1), (0, 1), (—1, 1)} (we omit the neighbors that 
can be obtained by reflection around (0,0)) and the interaction tensor as 



ki ii + i2 is even 
^2 h + ^2 is odd. 



Such material is illustrated in Fig. [2l 

This example was motivated by the study of Friesecke and Theil [20] , where a similar model was 
considered. Friesecke and Theil considered the model with springs similar to the one illustrated in 
Fig. O which however was nonlinear due to nonzero equilibrium distances of the springs (so that 
the energy of the spring between masses Xi and xj is proportional to \xi — Xj\'^ — Iq, where Iq is the 
equilibrium distance). They found that with certain values of parameters the lattice looses stability 
to non-Cauchy-Born disturbances and the lattice period doubles (thus the lattice ceases to be a 
Bravais lattice). 

The results, given with no details of actual derivation, are the following: The period of spatial 
oscillations in this case is (2, 2). The function x has the form x = xO^j) = ( — 1 )-?i +J2 ^^^^^^ ^ / (here 
/ is the 2x2 identity matrix). The homogenized tensors have the form 

= v4 = (^. + + + «^3) /. = = 

7.4 HQC 

In this subsection we sketch a formulation the HQC method based on the discrete triangular 
elements. Namely, we choose a partition T with triangles S 7~ of the original atomistic domain 
X = {i : (1, 1) < i < N}. The space of piecewise affine deformations is denoted as U, 



H 

per • 



Inside each triangle Sk choose a sampling rectangle S"^/^^ C Sk of the size pi x p2 atoms. Define 
the HQC energy variation: 

R 

Sk&T r=l 

where \Sk\ is the area of the triangle, and the microfunction IZk [w^) is defined on 5^°'' so that 
TZk {w^)-u^ G U^ieZ"^) and 

R 

J2 {ADrTZk {w^) , Drs)^^^grcp = Vs G [/^(eZ^). 

r=l 

As before, the function TZf^ (t/;^) can be extended on the whole triangle Sk if required. 
The variational problem to be solved thus becomes 



where 
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Figure 4: Strain Dui of the solution of the ID hnear problem: the schematically shown complete 
solution (left) and the closeup of the micro-structure for 31 atoms (right). 

8 Numerical Examples 

We solve numerically several model problems to illustrate the performance of HQC. We consider 
the linear and the nonlinear ID model problems (Sections 18.11 and I8.2p . followed by the two 2D 
linear model problems (Sections 18.41 and 18. 4p . We also study dependence of the numerical error on 
p, the spatial period of heterogeneity of the material (Section l8.3p . 

The aim of the numerical experiments is twofold. First, we verify numerically the sharpness of 
the obtained error for the ID linear case. Second, we investigate whether the HQC convergence 
results obtained for the nearest neighbor linear interaction in ID are valid for more general cases. 
The numerical results show that all theoretical conclusions made in Section [6] are also valid for 
finite range nonlinear interaction or for 2D problems. 

8.1 ID Linear 

In the first numerical example we solve the problem ([9j) for the linear interaction case with the 
period of spatial oscillation p = 2 and number of interacting neighbors R = 3. The potential is 
defined as ^ 

ipi,i+r{z) = -ki^i+r3^~''{z -rf (1 < r < R), 

where 

, f 1 i \s even 

[2 z IS odd 

Such potential is periodic, hence, as suggested by Theorem El Cmod = 0. The number of atoms was 
N = 2^^ = 16384, and the external force was taken as 

/i = sin(l + 2^Xi). 
The strain Dui for such problem is shown in Fig. 21 

Figure [5l is aimed to illustrate theorems HI and [5l It can be seen that the homogenized HQC 
solution converges to the exact solution in the L^-norm, does not converge in the ff^-norm, but the 
post-processed HQC solution does converge in the H^-noiia. The convergence of the homogenized 
solution in the L^-norm is exactly as suggested by Theorem [H first it decreases with the second 
order as H is refined, and later it stays constant as H is refined further. The 0(e) behavior of the 
lower bound of \\u^ — «||^2 is illustrated in Fig. [H 
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Figure 5: Results for the ID linear problem: errors of the post-processed HQC solution it and 
the homogenized solution in different norms. The errors behave in accordance with Theorems 
a and El 



error 
1 

0.1 - m 

0.01 - m □ N=2 



,n n n n □ D 



10"* 0.001 0.01 0.1 



H 



14 



0.001 - ^ + N=2'- 

o o o & 

10"' - + + + + g O N=2'« 



Figure 6: Results for the ID linear problem: errors ||u — n||^2 for N = 2 , N = 2 , and N = 2^. 
We can see that the plateau for small H follows the 0(e) = 0(A^^^) behavior, as predicted by 
Theorem HI 
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Figure 7: Results for the ID linear problem: errors of the post-processed HQC solution u '"^ and 
solution by the naive QC method u*^^. Here u is the exact solution, Av is the averaging operator 
defined as Av(u)i = !h+^i+l_ xhe graph illustrates that the naive of QC to a complex material 
fails, while the HQC successfully convergence to the exact atomistic solution. 



The errors of the post-processed HQC solution u ''^ and the solution by the naive QC method 
are shown in fig. [71 It can be seen that only the solution by HQC converges to the exact 
solution u, but the solution with the naive QC method does not converge, even when compared to 
the averaged exact solution (averaging operator is Av{u)i = ) or computed in an L^-norm. 

These findings are similar to calculations of Tadmor et al [33] which show that assuming a linear 
interpolation for a silicon crystal greatly overestimates its strain energy density. 



8.2 ID Nonlinear 

We solve the problem ([9|) for a general nonlinear interaction case, with the period of spatial oscil- 
lation p = 2, number of interacting neighbors i? = 3, and number of atoms = 2^^ = 16384. We 
chose Lennard-Jones potential 

/ \ -6 / \ -12 

^i,^+r{z) = -2 + 7^ {l<r<R) 

\H,i+r J \H,i+r / 



with the varying equilibrium distance 



1 z is even 

9/8 i is odd. 



The external force was taken as 

fi = 50sin(l -F27rXi) . 
The strain Dui for such problem is shown in Fig. [8j 

Figure [9| plots the errors of HQC in H^- and L^-norms. The results observed are qualitatively 
the same as the results on Fig. [7] obtained for the linear case. Also, the results of a naive application 
QC method to the nonlinear problem will be the same as shown on Fig. [7| for the linear case. Thus, 
we conclude that the convergence estimates derived for the linear case are also valid for the nonlinear 
case, as the conducted numerical experiments show. 
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strain 




Figure 8: Strain Dui of the solution of the ID nonhnear problem: the schematically shown complete 
solution (left) and the closeup of the micro-structure for 31 atoms (right). 
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Figure 9: Results for the ID nonlinear problem: errors of the post-processed HQC solution u '"^ and 
the homogenized solution in different norms. The errors behave in accordance with Theorems 
a and El 
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p = 2 p = 4 p = 8 p = 16 


1st test case (linear) 
2nd test case (nonlinear) 


0.040 0.043 0.044 0.042 
0.019 0.018 0.015 0.017 



Table 1: Dependence of the bound Cg = max-^^^ — (cf. Theorem [5]) on p. It can be seen that 
Cg essentially does not depend on p. 

8.3 Convergence for Different Periods p 

In our analysis of the equations and the computational method, we kept the dependence on p 
implicit, because derivation of estimates which are sharp w.r.t. p is much more technical. In this 
section we numerically address this issue. 

The first test case is similar to the one in Section 18.11 We fixed the bounds for /cj^j+r in (j84p 
between 1 and 2 and randomly generated the values of fcj^j+r with the periods p = 2, 4, 8, 16. We 
estimate the constant Cg in Theorem [S] as 

Cg = m ax , 

where the maximum is taken for H = 2^^, 2~^, . . . , 2""^. The second test case is similar to the one 
in Section 18.21 with h^i+r randomly generated between 1 and . 

The results for both test cases are shown in Table [H It can be seen that the constant Cg 
essentially does not depend on p. 

This finding is important in applications, for instance, to shape memory alloys that may change 
their crystalline structure in the course of loading/unloading. Motivated by such applications, 
the authors of [13] designed the adaptive strategy of choosing p, called Cascading Cauchy-Born 
kinematics, for the complex lattice QC method. They also presented an example of application 
of their method to the ID model problem exhibiting period-doubling bifurcations. The present 
findings indicate that increase of period p does not aff'ect the accuracy of the method. 

Independence of error bounds on p is also important for modeling amorphous materials, such 
as glasses or polymers. Amorphous materials do not have a spatial period, instead they exhibit 
some random structure. By analogy with application of numerical homogenization to PDEs with 
random tensors, one could take p large enough to capture the variation of the microscopic structure 
of the amorphous material, and expect that it will not affect the accuracy of representation of the 
macroscopic deformation as the mesh size H is refined. 



8.4 2D Test Case 1 

We consider the example of material discussed in Subsection 17.3.31 with e = 2^^^, A''i = A''2 = 2^^, 
ki = 1, /c2 = 2, A;3 = 0.25, 

f iri — cos(7rji/A''i)^— cos(7rj2 

/7V2)2 / sin(27rfi/iVi) \ _ j 
~ V sin(27ri2/iV2) ) ^' 

where / is determined so that the average of fi is zero. The total number of degrees of freedom 
of such system is approximately 8 • 10^. The solution for such test case is shown in fig. [10] (the 
illustration is for A^i = A'^2 = 64). 
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Figure 11: Illustration of a 2D triangulation. Larger atoms comprise sampling domains for HQC. 
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Figure 12: Results for the 2D test case 1: error depending on the mesh size H. The L^-error of 
the homogenized solution and the i?^-error of the post-processed solution u^''^ are shown. The 
errors behave in accordance with the ID analysis (Theorems [4] and [5]) . 



The atomistic domain is triangulated using nodes and K = 2t^ triangles (t = 2,4, . . . ,2^"). 
In each triangle a samphng domain is chosen, each samphng domain contains four atoms (see 
ihustration in fig. [TTI) . The number of degrees of freedom of the discretized problem is 2t^ . 

The error of the solution for different mesh size H (H = 0.5, 0.25, . . . , 2"^^^) is shown in fig. [T2j 
The results are essentially the same as in ID case: the method convergences with the first order 
of mesh size in the i7^-norm and with the second order in the L^-norm. We also see the plateau 
for the L^-error of the homogenized solution u'^'^ . It is remarkable that all the conclusions of ID 
analysis (cf. Theorems [H and [5]) are also valid for the 2D computations. 



8.5 2D Test Case 2 

The second test case is analogous to the previous one, but with the different tensors tpr describing 
the atomistic bonds. The tensors ipr were chosen to have the following (randomly generated) values: 



V'(l,0),i 



1.3 




even, 12 even 


1.6 


k 


even, ^2 odd 


1.8 


k 


odd, i2 even 


1.2 


k 


odd, 12 odd. 


0.3 


k 


even, 12 even 


0.8 


k 


even, ^2 odd 


0.6 


k 


odd, i2 even 


0.4 


k 


odd, 12 odd. 



^(0,1),' 



1.5 


k 


even. 


12 even 


1.7 


k 


even. 
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For a tensor with such a random structure, the homogenized tensor can only be precomputed 
numerically, and in the case of a nonlinear problem should be found in the course of the actual 
computation. The solution (for A^^i = N2 = 64) is shown in fig. [T3j 

The error of the solution for different number of degrees of freedom is shown in fig. [TH The re- 
sults are similar to the results of all the previous test problems. Again, the results are in accordance 
with Theorems |4] and [5l 
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Figure 13: Atomic equilibrium configuration for A'^i = = 64 for the 2D test case 2. Deformation 
of the whole material (left) and a close-up (right). 
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Figure 14: Results for the 2D test case 2: error depending on the mesh size H. The L^-error of 
the homogenized solution and the i?^-error of the post-processed solution u^'^ are shown. The 
errors behave in accordance with the ID analysis (Theorems H] and [5]) . 
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9 Summary and Conclusion 

We have considered the problem of modehng materials with complex atomistic lattice. We have 
proposed a discrete homogenization framework to analyze the QC method for complex crystalline 
materials. This framework allowed us to prove convergence (in ID) for the QC method proposed in 
[33]. Numerical homogenization has also been used to formulate the QC method. The equivalence 
of this algorithm to the QC method of [3_3j is discussed in detail in [3]. We have also shown how 
to apply the presented technique in a 2D setting. The ID and 2D numerical examples presented 
verify validity of the analysis in more general setting. We note that the extension of the algorithm 
proposed in this paper to simulate atomistic materials at finite temperature or non-crystalline 
materials is of high interest. This is a topic for future research. 
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Appendix A 

Lemma 10 (Discrete Poincare inequality 1). Let g € and X^^i = 0. Then 



L ^2 L-l 



6 

i=l 1=1 

Proof. We start with noticing that lemma A.l in [27, p. 87] applies to g and states that 

L-1 



where 



- j <i 



Then with the help of the Cauchy-Schwarz inequality one obtains 

E < E f E i^.+i - ^.i^^..) ^ E f E i^.+i - f E ^l- 1 ' 
i=i i=i \j=i J i=i \j=i J \j=i 

where by direct computation 

L L-1 L-1 L L-1 ( , ..2 / ■\2 \ 

EE<,= EE*L = E (^) i+(i) (i-i) 

i=i j=i j=i i=i j=i \ ^ ^ ^ / 

^ ^^ (L-j)j _L^-l 

^ L 6 6 ■ 



□ 
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Lemma 11. Let g £ (IR^)^.- Then 



'# 

L 

2L 



^L + l-2k . . 
^ 2^ of \9i-k+l - Qi-k)- 



k=l 

Proof. Direct computation of the right-hand-side (RHS) yields: 

L 

2L 



RHS — ^ — ^t^jy {gi-k+i - Qi-k) 

k=l 



L + 1 - 2(/c - 1) A L + 1 - 2A; 

fc=o fc=i 
L - 1 4-^2 L - 1 



'9i - 9i-k H — :Tj—9i-L- 



2L ^ 2L " 2L 

fc=i 

L-l 

Notice that due to periodicity gi^i = gi and due to average of gf being zero, Yl 9i-k = —9i- Hence 

fc=i 

T.TTO 2 L-l 

□ 

If we consider the periodic extension of the sequence then the estimate of lemma [10] will have 
a slightly better constant: 

Lemma 12 (Discrete Poincare inequality 2). Let g G Then 

^\9i? < -^^\9i+i - 9i? ■ 

i=l i=l 

Proof. Then by using lemma [TT] and the Cauchy-Schwarz inequality one obtains 

V-i 12/ ^ (^L + l-2k. X 

2^\9i\ < 2^ 2^ ^ {9i+i~k - 9i-k) \ 

i=l i=l \k=l / 



\2 

-1-k — Qi-k) 



i=l k=l ^ ' k=l 

- 1 

^i9i+i-k - Qi-kf 

i=l k=l 

L"^ - 1 ^ ^ 
—[^ Y.^9j+i - 9j? < Y.^9j+i - 9jf- 



□ 
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Corollary 1 (Discrete Poincare inequality for U^). The functional | • (cf. ^) defines a norm 
on . For u G the following inequality holds: 



Ml2 < :r^mm- 



Lemma 13 (Inverse discrete Poincare inequality). For u G [/^"p the following inequality holds: 

e\u\y^n,q < 2\\u\\li (86) 

for 1 < q < oo. 
Proof. 

□ 

Lemma 14. For u G the following inequality holds: 

Proof. Using the discrete Poincare inequality yields: 

{u,v)- 1 

\u\^~l = sup -j — j < sup -= = p||u||i2. 



□ 



